Designing reserve networks requires considering both ecological and socioeconomic factors. Marxan is a commonly used spatial planning tool that quantifies these factors to generate efficient reserve networks. Here, we will use species and parcel data from a 2014-15 Bren Group Project in the Morro Bay Watershed to examine the optimate regional reserve network.
Prioritizr is an R package that, like Marxan, can guide systematic reserve design by solving conservation prioritization problems. In this tutorial, we will be using the marxan_problem() function which can read Marxan input data. More information on the prioritizr package can be found at https://mran.microsoft.com/snapshot/2018-03-04/web/packages/prioritizr/vignettes/prioritizr_basics.html.
A really useful tutorial in using prioritizr, specifically in using summed solutions, can be found here: https://cran.r-project.org/web/packages/prioritizr/vignettes/tasmania.html. This resource helped frame this tutorial.
For the solve() function, you must have a solver installed. For the the process with summing and multiple runs process, this must be gurobi. This requires following the instructions below.
install.packages('/Library/gurobi902/mac64/R/gurobi_9.0-2_R_3.6.1.tgz', repos=NULL) (if on a mac and downloaded the most recent version of gurobi, don’t download newest version of R! (4.0.0; at least as of 5/2020 this was not functional))We will be using a case study (and the data they prepared) from a past group project on the Morro Bay National Estuary Program.
Species data
Parcel data
library(here)
library(janitor)
library(prioritizr) # marxan
library(sf) # spatial features
library(slam) # to use gurobi
library(gurobi) # solver
library(ggmap) # basemaps
library(patchwork)
library(tidyverse)
If you do not have any of these packages installed, run the following in the console: install.packages("package name")
Note: this will not work for installing gurobi. Instructions are in the introduction section.
All .csv files in this tutorial are unaltered from the .xlsx files provided in the R drive. They have only been renamed and saved as csvs. The “data” folder is identical to the MorroBay_data folder provided in the R drive, it has only been renamed.
spec <- read_csv(here("data", "morro-bay-spec.csv")) %>%
head(140) %>% # mine reads in extra blank rows - skip if yours does not
rename("amount" = "target") # target is called "prop" (relative) or "amount" (absolute)
pu <- read_csv(here("data", "morro-bay-pu.csv")) %>%
select(1:3) # mine reads in an extra blank column - skip if yours does not
puvsp <- read_csv(here("data", "morro-bay-puvspr.csv")) %>%
select(1:3) %>%
head(11849)
status <- read_csv(here("data", "spec-name-status.csv")) %>%
select(1:3) %>%
head(140)
parcels <- read_sf(dsn = here("data"), layer = "MorroBay_parcels") %>%
clean_names()
ggplot(data = parcels) +
geom_sf()
The function marxan_problem() in the prioritizr package gives us the “canned” approach that works for our purposes. If more customizations are desired, feel free to explore the problem() function instead.
Arguments for marxan_problem (for more information, see ?marxan_problem):
x = pu)spec = spec)puvspr = puvsp)bound = NULL)Also note that in solver, we will use 100 runs (number_solutions = 100)
marxan_1 <- marxan_problem(x = pu,
spec = spec,
puvspr = puvsp,
bound = NULL)
marxan_1_problem <- marxan_1 %>%
add_gurobi_solver(gap = 0.15) %>%
add_pool_portfolio(method = 2, number_solutions = 100)
# gap = optimality gap. 0.15 gives us results closest to that of using marxan software. Learn more under ?add_gurobi_solver
# method 2 finds a specified number of solutions that are nearest to optimality. Learn more under ?add_pool_portfolio
marxan_1_soln <- solve(marxan_1_problem)
## Gurobi Optimizer version 9.0.2 build v9.0.2rc0 (mac64)
## Optimize a model with 140 rows, 670 columns and 11849 nonzeros
## Model fingerprint: 0xd2c5bccb
## Variable types: 0 continuous, 670 integer (670 binary)
## Coefficient statistics:
## Matrix range [1e+00, 1e+00]
## Objective range [2e+00, 5e+06]
## Bounds range [1e+00, 1e+00]
## RHS range [3e-01, 2e+02]
## Found heuristic solution: objective 4.238988e+07
## Presolve removed 116 rows and 205 columns
## Presolve time: 0.01s
## Presolved: 24 rows, 465 columns, 3066 nonzeros
## Variable types: 0 continuous, 465 integer (465 binary)
## Presolve removed 1 rows and 0 columns
## Presolved: 23 rows, 465 columns, 3058 nonzeros
##
##
## Root relaxation: objective 2.840941e+07, 14 iterations, 0.00 seconds
##
## Nodes | Current Node | Objective Bounds | Work
## Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
##
## * 0 0 0 2.840941e+07 2.8409e+07 0.00% - 0s
##
## Optimal solution found at node 0 - now completing solution pool...
##
## Nodes | Current Node | Pool Obj. Bounds | Work
## | | Worst |
## Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
##
## 0 0 - 0 - 2.8409e+07 - - 0s
## 0 0 - 0 - 2.8409e+07 - - 0s
## 0 2 - 0 - 2.8409e+07 - - 0s
##
## Explored 1958 nodes (383 simplex iterations) in 0.15 seconds
## Thread count was 1 (of 4 available processors)
##
## Solution count 100: 2.84094e+07 2.84097e+07 2.842e+07 ... 3.34128e+07
##
## Optimal solution found (tolerance 1.50e-01)
## Best objective 2.840941400000e+07, best bound 2.840941400000e+07, gap 0.0000%
marxan_1_ssoln <- marxan_1_soln %>%
mutate(sum = rowSums(.[6:105])) %>%
select(id, cost, status, locked_in, locked_out, sum)
hist(marxan_1_ssoln$sum,
main = "Selection Frequencies",
xlab = "Number of runs that the unit was selected")
parcels_marxan_1 <- inner_join(parcels, marxan_1_ssoln, by = "id")
ggplot(data = parcels_marxan_1) +
geom_sf(aes(fill = sum),
color = "white",
size = 0.05) +
scale_fill_gradient(low = "slategray2",
high = "navy") +
labs(fill = "Summed \nSolution") +
theme_minimal()
This step is optional, but will make your map look better and is good to learn for further spatial analyses in R. Unfortunately, the options for basemaps in R are limited. Here, I have used the ggmaps package because it is the most accessible, however, it is still not very accessible compared to other R packages because you’ll need to get an API key. Instructions are here: https://cran.r-project.org/web/packages/ggmap/readme/README.html
For a ggmap cheatsheet, including the different basemaps in ggmaps, see: https://www.nceas.ucsb.edu/sites/default/files/2020-04/ggmapCheatsheet.pdf
morrobay <- get_map(location = c(lon = -120.7665, lat = 35.335),
zoom = 12,
maptype = "terrain-background", # background means no references, omit if want references
source = "google")
all_spec <-
ggmap(morrobay) +
geom_sf(data = parcels_marxan_1,
aes(fill = sum),
color = "white",
size = 0.1,
alpha = 0.85,
inherit.aes = FALSE) +
coord_sf(crs = st_crs(4326)) +
scale_fill_gradient(low = "slategray2",
high = "navy") +
labs(title = "All Species",
fill = "Summed \nSolution",
x = NULL,
y = NULL) +
theme_minimal()
all_spec
Here we create spec and puvsp data frames similar to those we had for all species, but containing only endangered or threatened species.
end_status <- status %>%
mutate("endangered" = case_when(
str_detect(status, pattern = "endangered") == TRUE ~ "yes",
str_detect(status, pattern = "threatened") == TRUE ~ "yes",
T ~ "no")) %>%
filter(endangered == "yes")
end_spec <- merge(end_status, spec, by = "id") %>%
select(id, amount, spf, name.x) %>%
rename("name" = "name.x")
puvsp_id <- puvsp %>%
rename("id" = "species")
end_puvsp <- merge(end_status, puvsp_id, by = "id") %>%
rename("species" = "id")
marxan_end <- marxan_problem(x = pu,
spec = end_spec,
puvspr = end_puvsp,
bound = NULL)
marxan_end_problem <- marxan_end %>%
add_gurobi_solver(gap = 0.15) %>%
add_pool_portfolio(method = 2, number_solutions = 100)
marxan_end_soln <- solve(marxan_end_problem)
## Gurobi Optimizer version 9.0.2 build v9.0.2rc0 (mac64)
## Optimize a model with 18 rows, 670 columns and 970 nonzeros
## Model fingerprint: 0x6e88f508
## Variable types: 0 continuous, 670 integer (670 binary)
## Coefficient statistics:
## Matrix range [1e+00, 1e+00]
## Objective range [2e+00, 5e+06]
## Bounds range [1e+00, 1e+00]
## RHS range [9e-01, 1e+02]
## Found heuristic solution: objective 4.174920e+07
## Presolve removed 14 rows and 203 columns
## Presolve time: 0.00s
## Presolved: 4 rows, 467 columns, 495 nonzeros
## Variable types: 0 continuous, 467 integer (467 binary)
## Found heuristic solution: objective 3.534249e+07
## Presolved: 4 rows, 467 columns, 495 nonzeros
##
##
## Root relaxation: objective 2.677119e+07, 1 iterations, 0.00 seconds
##
## Nodes | Current Node | Objective Bounds | Work
## Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
##
## * 0 0 0 2.677119e+07 2.6771e+07 0.00% - 0s
##
## Optimal solution found at node 0 - now completing solution pool...
##
## Nodes | Current Node | Pool Obj. Bounds | Work
## | | Worst |
## Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
##
## 0 0 - 0 - 2.6771e+07 - - 0s
## 0 0 - 0 - 2.6771e+07 - - 0s
## 0 2 - 0 - 2.6771e+07 - - 0s
##
## Explored 4914 nodes (1716 simplex iterations) in 0.41 seconds
## Thread count was 1 (of 4 available processors)
##
## Solution count 100: 2.67712e+07 2.67732e+07 2.67736e+07 ... 3.14836e+07
##
## Optimal solution found (tolerance 1.50e-01)
## Best objective 2.677119300000e+07, best bound 2.677119300000e+07, gap 0.0000%
marxan_end_ssoln <- marxan_end_soln %>%
mutate(sum = rowSums(.[6:105])) %>%
select(id, cost, status, locked_in, locked_out, sum)
hist(marxan_end_ssoln$sum,
main = "Selection frequencies",
xlab = "Number of runs that the units were selected")
parcels_marxan_end <- inner_join(parcels, marxan_end_ssoln, by = "id")
ggplot(data = parcels_marxan_end) +
geom_sf(aes(fill = sum),
color = "white",
size = 0.05) +
scale_fill_gradient(low = "slategray2",
high = "navy") +
labs(fill = "Summed \nSolution") +
theme_minimal()
end_spec <-
ggmap(morrobay) +
geom_sf(data = parcels_marxan_end,
aes(fill = sum),
color = "white",
size = 0.1,
alpha = 0.85,
inherit.aes = FALSE) +
coord_sf(crs = st_crs(4326)) +
scale_fill_gradient(low = "slategray2",
high = "navy") +
labs(title = "Endangered & Threatened Species",
fill = "Summed \nSolution",
x = NULL,
y = NULL) +
theme_minimal()
end_spec
This step isn’t necessary, but is really convenient if you want to present your maps together in your report. patchwork uses PEMDAS to combine ggplots together into one graphic. For example, plot_1 + plot_2 will make an image of the plots side by side, whereas plot_1 / plot_2 will make an image of plot_1 over plot_2. For more information, see https://github.com/thomasp85/patchwork.
First, make the all_spec graph without a legend so that the combined image has only one legend.
all_no_lgnd <- all_spec +
theme(legend.position = "none")
Next, combine the graphs with patchwork
spec_graphs <- all_no_lgnd + end_spec
spec_graphs
Finally, save the image to add to your report!
ggsave("spec-graphs.png")